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ABSTRACT 


A  method  is  described  in  which  difference  equations  are  derived  by- 
application  of  the  principle  of  virtual  work  to  a  particular  fluid  model. 
Some  other  innovations  involve  the  viscous  pressure,  different  time 
intervals  for  different  points,  and  a  method  of  handling  the  collision 


of  free  interfaces 


These  notes  describe  several  innovations  in  Lagrangian,  time- 
dependent,  two-dimensional  hydrodynamics.  These  are  incorporated  in  a 
code  called  WAT  now  in  use  on  the  IASL  IBM  Model  70^  computers .  When 
applied  to  problems  that  involve  only  one -dimensional  motion  it  appears 
to  be  as  accurate  as  one-dimensional  treatments  with  comparable  mesh 
size.  In  several  two-dimensional  problems  it  has  given  reasonable  results. 

I .  The  Fluid  Model  —  Differential  and  Difference  Equations 

We  partition  the  fluid  into  a  number  of  discrete  zones,  adopt  a 
specific  physical  model  of  the  partitioned  fluid,  and  use  D'Alembert’s 
principle  to  obtain  exact  equations  of  motion  for  the  model.  These  are 
differential  equations  in  time,  and  we  difference  them  in  a  straight¬ 
forward  way. 

The  fluid  is  divided  into  cells  of  quadrilinear  cross  section,  as 
in  Figure  1.  We  will  take  all  inertia  to  be  concentrated  on  point  par¬ 
ticles  at  the  cell  vertices;  the  material  within  each  cell  exerts 
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pressure  determined  through  an  equation  of  state  by  its  internal  energy 
E  and  volume  V. 

Suppose  point  (i,j)  is  displaced  infinitesimally  by  S_s.  Work  is 
done  on  the  surrounding  four  cells  and  the  system's  energy  is  modified 
by  6E. 


6E  =  PlJ  8s  •  Wij+Pij.jSS-Wij.1  + 


+  p.  .  .5s  »W.  . . 

1-1,  J - i-lj 


(1) 


According  to  a  general  principle  of  mechanics  (D'Alembert's 
principle)  a  force  F  acts  on  the  particle  (i,j) 


5E  =  -5£ 


or 


F.  .  =  p.  .W.  .+p.  .  -W.  .  ,+p.  ,  .  ,W.  ,  .  , 


(m 


The  equation  of  motion  of  point  (i,j)  is  then 


r 

£iJ  “  Mij 


(2) 
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which  may  be  differenced  in  time 


t+At  0  t  t-At 

£i.i  + 

(At)2 


F*. 
_  ZM 
•  ”  M.  . ' 
ij 


(3) 


In  addition  we  need  to  keep  account  of  the  internal  energy  per  unit 
mass  through 


de 


p.  .dV.  . 

a  ^  3-J 

^  (PDVo) 


(M 


where  (p  V  )  is  the  mass  of  fluid  in  the  cell, 
o  o 

We  take  the  mass  M. .  of  the  particle  at  vertex  (i, J )  to  be  half  the 
3*  J 

mass  of  the  shaded  regions  in  Figure  1. 

The  difference  equations  are  written  in  full  in  the  appendix  for 
cylindrical  (r,z)  coordinates.  It  is  a  straightforward  task  to  compute 
them  for  other  coordinate  systems.  In  cartesian  (x,y)  coordinates  they 
reduce  to  Kolsky's  equations2  only  for  the  case  that  the  mesh  is  exactly 
rectangular . 

Kolsky’s  equations  are  obtained  by  regarding  the  vertices  as  re¬ 
presentative  points  fixed  in  a  continuous  fluid.  A  pressure  gradient 
at  the  vertex  is  estimated  and  used  to  compute  an  acceleration.  This 
can  be  done  in  many  ways,  and  the  present  method  can  be  regarded  as  just 
another,  different,  way  of  doing  this .  However,  I  would  suggest  that  in 
a  highly  distorted  mesh,  where  it  is  difficult  to  define  a  sensible 
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pressure  gradient,  the  present  method  still  makes  physical  sense  --  the 
force  on  each  particle  is  determined  by  the  reaction  of  the  mesh  to  its 
displacement . 

II.  Artificial  Viscosity 

To  spread  shock  fronts  over  several  zones,  ve  introduce  an  artificial 
viscous  pressure  q,  following  Richtmyer  and  von  Neumann'*'.  A  somewhat 
detailed  interpretation  of  the  viscous  pressure’s  role  is  useful  in 
finding  a  form  appropriate  to  two  dimensions.  For  example,  our  exper¬ 
ience  is  that  the  isotropic  form,  q  proportional  to  (-r^r)2,  is  not 
adequate  in  a  mesh  far  from  square,  where  shocks  coming  from  different 
directions  see  cells  of  quite  different  thickness. 

In  one  dimension,  Richtmyer  and  von  Neumann’s  form  for  q  can  be 
written 

q  =  ap(Au)2  (5) 

where  p  is  the  density  in  the  cell,  a  is  a  dimensionless  constant  of  the 
order  of  unity,  and  Au  is  the  velocity  of  approach  of  the  two  boundaries 
of  the  cell. 

If  Au  is  a  good  deal  smaller  than  sound  velocity  in  the  cell,  no 
discontinuous  behavior  is  expected  and  numerical  integration  should 
proceed  smoothly.  Addition  of  q  to  the  pressure  can  be  regarded  as  an 
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expedient  to  increase  the  sound  velocity  in  the  cell  so  that  this  con¬ 
dition  is  always  met. 

To  see  this,  let  the  equation  of  state  of  the  fluid  be  represented 
by  a  7  law. 


P  =  (7-l)ep 

and  we  have  the  energy  equation 

de  =  (pfq)  ^|. 

P 

The  sound  velocity  is 

■  0.  -  0c  +  <&,<& 

a)f 

which  can  be  written,  using  Eq.  (5) 

c2  =  c2  +  a(7-l)(Au)2  (6) 

where  c  is  the  sound  velocity  with  q=0. 

2 

So,  one  sees,  the  effect  of  q  is  to  add  (Au)  multiplied  by  a 
factor  of  the  order  of  1  to  the  square  of  the  physical  sound  velocity. 


8 


and  this  will  insure  that  two  boundaries  of  a  cell  never  approach  each 
other  rapidly  compared  with  the  effective  sound  velocity  in  the  cell. 

Following  this  interpretation,  we  have  made  q  depend  on  the  largest 
of  the  (6)  velocities  of  approach  of  pairs  of  vertices  of  each  cell,  i.e.. 


q  =  a- 


KV 


~(Au«r)^~] 


some  Au*r  <  0 


max.  of  6  pairs 


(7) 


=  0  every  Au»r  >  0 

where  Au  and  r  represent  the  relative  coordinate  and  velocity  of  a  pair 
of  vertices. 


III.  Stability 

A  few  experiments  with  simple  stability  criteria  using  the  sound 
transit  time  through  an  average,  or  alternatively,  the  smallest,  dimen¬ 
sion  of  a  cell  have  convinced  me  of  the  value  of  a  more  accurate,  if 
more  complicated  criterion.  When  the  mesh  became  distorted,  only  a 
criterion  using  the  smallest  distance  seemed  adequate,  and  this  wasted 
a  great  deal  of  computation  by  nearly  always  demanding  too  short  an 
interval  —  in  the  limit  when  two  points  coalesce,  for  example,  a  zero 
interval  is  demanded. 

We  obtain  a  more  accurate  criterion  by,  as  usual,  linearizing  the 
equations  of  motion  about  an  average  motion,  following  von  Neumann  and 
Richtmyer^ .  Let  5z^,  Sr^  denote  a  (small)  departure  from  the  exact 
solution  of  the  difference  equation  (3).  The  pressure  in  each  cell  will 
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change  to  first  order  in  6z,  5r,  by 


5p  =  frc2  21  [S 5z  +  §? 6r]  (8) 

00  tn 

vertices 

2 

where  c  is  the  sound  velocity,  c  =  and  the  sum  egresses  the  total 
change  in  V  due  to  displacements  at  all  four  of  its  vertices.  For  the 
moment,  we  neglect  the  dependence  of  q,  the  viscous  pressure,  on  the  5’s. 
Then  we  obtain  equations  of  motion  for  the  6*s, 


.  t+At  „  t  _  t-At 
5zi.i  -  25Zi.i  +  6Zi.i 

At2 


dV,  ,  rdV 


M.  .vp  V  ~  ^ij 
ij  o  o  0  ij 


6z  . 
Zij  ij 


dV.  . 
+  ^ 


dV. 


^ 8riJ  +  dzi,ti  5Zi’J+1  +  ari’^,'1 


dV.  . 


avil  dVij 

8Zi+1.^1  +  dr!+i,jn  Sri+1.Jtl  +  3zi+i,j  i+1-J 


d'l.  , 

+  —  6r.Ll  ,  \  +  similar  terms  for  neighboring  cells 


1+1,  j 


:i+l,j) 


(1-1, j);  (i-l,j-l);  (i,j-l). 


(9) 


There  is  a  similar  equation  for  6r.  In  a  compressed  notation  for  the 
right  hand  side. 
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tf At  „  t  R  t-At 
5z. .  -  25z  +  5z 

_ 

(At)2 


=  D6t. 

z 


(9* ) 


An  accurate  stability  limit  would  be  obtained  by  finding  the  larg¬ 
est  eigenvalue  of  the  (linear)  operator  D 


DSz ..=  of  6z . 
max 


The  solutions  of  difference  equation  (9)  change  character,  from 
oscillatory  (neutrally  stable)  to  exponentially  growing  (unstable)  as  At 
is  increased  beyond  At^  , 


(At  )2  = 
v  max' 


(10) 


0) 


max 


We  construct  a  local  stability  limit  for  each  point  by  replacing  all 


bv  P  2  c)V  2  2 

of  the  factors  3—  pc,  pc  by  the  largest  of  them. 

dz  dr  dV  dV 

We  also  use  the  larger  of  the  two  factors  . 

1 0  1 0 


The  largest 

eigenvalue  will  then  occur  for  that  mode  in  which  all  terms  in  the 


approximate  D  are  in  phase,  and  we  construct  a  local  At^y  from  an  cd^ 


defined  by 


2 

0.  . 

u 


- 

^V. 

_ 

“  M,  . 
ij 

larger  ofS 

dz.  . 

► 

fev. . 

I 

Idr.  . 
^  ij- 

. 

(largest  dv, 

<at  any 
[vertex 


=  it-b 


'i  y 


(11) 


11 


Another  source  of  instability  is  the  viscous  pressure  q,  which  de¬ 
pends  on  a  relative  velocity,  Eq .  (7 ) .  To  first  order  in  the  perturbation 
a  term  will  be  added  to  the  right  hand  side  of  Eq.  (9)  for  each  cell 
having  (i,  j)  as  a  vertex.  The  one  added  to  the  terms  written  out  in 
Eq.  (9)  is 


.  dv.  .  a(p  V  ). . 
1  i,i  VHo  o'i.i 

m  . .  3z7T  v . . 

ij  ij  ij 


..t  .  t-At 

)  +  g|Au|^ 


AlSr'-Sr^*) 
At 


]• 


As  in  Section  II,  A  identifies  the  largest  rate  of  approach  of  any 


pair  of  points  around  the  cell.  Again,  we  let  all  coefficients  of  the 


6*s  take  on  the  value  of  the  largest  one,  and  assume  the  mode  that  makes 

2 

all  of  them  add.  The  result  is  that  we  add  to  o>^,  Eq.  (11): 


kb*1 

Zil.Ai 2 

At  "  M.  .  At 
IJ 


5H71 


larger  of< 


ij 

dVij 


ij-'  J 


largest  value 
among  the  cells 
of  which  (i,j) 
is  a  vertex 


(12) 


And  finally,  At^^  is  determined  for  each  point  by 


(At)2  b°  +  Atb^.  <  1. 
ij  IJ 

IV.  Use  of  Many  Time  Intervals 


(13) 


In  many  problems,  a  large  reduction  in  computation  can  be  realized 
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by  moving  each  point  vith  the  time  interval  demanded  for  its  stability. 

In  this  section  the  scheme  used  in  WAT  to  do  this  is  outlined. 

Up  to  sixteen  different  time  intervals  are  defined,  each  a  factor 
two  smaller  than  its  predecessor.  Each  point’s  acceleration  is  then 
advanced  with  the  largest  interval  that  meets  its  stability  criterion; 
its  acceleration  is  taken  constant  over  its  interval,  but  its  position 
is  advanced  as  needed  to  advance  the  pressure  in  cells  that  touch  a  point 
being  advanced  with  a  smaller  time  interval.  The  velocity  is  always 
advanced  to  the  center  of  the  interval  through  which  the  position  is  to 

O 

be  advanced;  positions  are  then  always  accurate  to  order  (At)  .  A  table 
(t?.)  is  kept  current  of  the  time  to  which  each  position  has  been  advanced, 

X  J 

and  another  (tT .)  of  the  state  of  the  velocities.  In  Figure  2  an  outline 

X  J 

flow  chart  of  the  calculation  is  shown.  The  time  intervals  are  At^, 

=  0,1,2,  ;  At  is  the  largest  interval.  F  and  F*\  are  flags 

max  o  co  x  j 

(up  or  down)  that  tell,  respectively  whether  a  given  interval  has  just 
been  processed,  and  whether  a  given  cell's  pressure  has  been  brought  up 
to  the  current  time.  A  table  is  available  of  the  i,j  coordinates  of  all 
points  to  be  moved  with  a  given  interval  At^. 

V.  Collision  of  Free  Surfaces 

WAT  is  provided  with  a  routine  for  handling  the  collision  of  two 
free  surfaces.  In  figure  3  a  typical  situation  is  depicted,  where  the 
points  and  cells  have  been  numbered  one  dimensionally. 
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FjTj/  up? 


no 


Advance  coordinates  to  t 
compute,  volume  advance, 
p,e 


Figure  2 


Ik 


Figure  3 

« 
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The  first  problem  is  to  compute  the  time  at  which  a  point  on  one 
free  surface  crosses  the  other,  as  point  j'  has  done.  We  make  use  of 
the  sense  of  a  point  with  respect  to  each  line  segment  defining  the  sur¬ 
face  that  it  may  cross;  it  defines  a  point  as  being  "inside"  or  "outside" 
of  the  region  bounded  by  the  line.  The  sense  of  point  j  with  respect  to 
the  line  through  j 1 ,  j'+l  is  the  sign  of  the  expression 

(w)(rj'+i-V  -  (rrrj')(zo’+i'zj')- 

A  point  has  not  crossed  a  given  segment  unless  its  coordinates  lie  between 
those  of  the  end  points  of  the  segment,  within  a  tolerance  given  by  the 
travel  of  the  point  during  a  time  cycle,  and  unless  its  sense  with  respect 
to  the  segment  has  changed.  When  both  of  these  conditions  are  fulfilled, 
we  compute  the  time  and  coordinate  of  the  intersection  of  the  path  of  the 
particle  with  the  segment.  In  general,  one  gets  two  solutions  for  each 
segment,  and  further,  the  above  conditions  may  have  been  met  for  more 
than  one  segment,  so  the  intersection  is  rejected  unless  the  time  of 
intersection  is  less  than  At,  and  unless  the  sense  of  the  path  —  from 
"outside  to  inside"  —  is  correct  at  the  intersection. 

In  general  it  is  necessary,  for  each  point  j,  to  examine  every  seg¬ 
ment  j',  j’+l  of  the  opposing  interface.  The  labor  is  reduced  in  WAT  by 
first  locating  all  of  the  points  in  a  coarse  Eulerian  mesh.  Then  only 
if  two  points,  not  neighbors,  fall  within  the  same  or  neighboring  cells 
does  the  code  look  for  intersections.  Possible  crossings  by  each  point 
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of  every  segment  connected  to  one  of  the  points  are  examined. 

When  a  valid  intersection  is  found,  the  three  particles  involved  — 

the  intersecting  one  and  the  two  that  define  the  line  segment  —  are 

taken  to  collide  inelastically.  That  is,  the  center  of  mass  velocity 

is  given  to  each  of  them,  and  the  relative  kinetic  energy  is  distributed 

* 

as  internal  energy  among  the  participating  cells . 

From  then  on,  point  j  is  attached  to  the  interface  of  cell  j  •;  that 
is,  the  volume  of  cell  j‘  is  computed  using  the  segments  j*,j  and  j,j’+l 
as  its  boundary,  and  the  equation  of  motion  of  point  j  has  a  term 

V)^S’ 

added  to  the  force.  It  is  necessary  to  maintain  two  tables:  for  each 
particle  j,  the  cell  to  which  it  is  attached,  if  any;  for  each  cell  j • , 
the  foreign  point  or  points  that  are  along  its  inner  boundary,  if  any. 


Other  models  for  the  inelastic  collision  were  tried,  but  this  one  led 
to  the  smoothest  subsequent  motion  of  the  collision  interface. 
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APPENDIX 


Difference  Equations  in  Cylindrical  Coordinates 

The  system  is  assumed  rotationally  symmetric  about  the  r=0  axis. 
Indexing  is  as  in  Figure  1. 

volume  Vi(j  = 

"  (zi+l,j‘Zij)(ri,>rrij)] 

+  f(ri+lj4-l+ri+lj+ri,H-l^  [(zi+lj"zi+l^l^ri^l"ri+l«^l^ 
"  ^ziJH"zi+lJH^ri+lj“Pi+l^l|l 

acceleration 
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(r nfl-2rn  +rn_1) 

iu  (  SoiJ 


ir  ij 


=  P 


At 


i-l,j 


(ri j+ri-l, j+ri, (zi, j+l"Zi-l, 0^ 
+  *Zi-l,  j"zij^ri,Jtl"rij* 


+  P. 


ij 


(rij+ri,jfl+ri+l,j)(zi+l,j"Zi,>l) 
+  ^zi,jfrzij^Pi+i,3"Pij^ 

[■  (zi+l,j'zij)(ri,j4-l“rij^ 


+  p. 


i,j-l 


trij+ri+l,j+ri,j-l^zi,j-rzi+l,j) 
+  (zi+l,j“zij)(ri,j-l-rij) 

“  (zi,j-l"Zij)(ri+l,o”rij) 


(ri j+ri, J-1+Pi-1, j } (zi-l, j“Zi, J-l ) 
+  <zi,j-i“zit5)(ri-i,rrij) 


-  (zi-i,j-zij)(ri,j-i-rij) 
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.  trf  At  _  t  ,  t-Ati 

1  M  (zij  “2zij  Zij  _  Tr  +r  +r  1  lr  -r  1 


7T  ij 


At 


pij  [rij+ri,>l+ri+l,j]  [ri,^l‘ri+l,J 
Pi,j-1  [rij+ri+l,j+ri,  j-l]  [ri+l,j“ri,J-l] 
Pi-1,  J-l  [rU+ri,J-l+ri-l, i] 
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